function [K_mod_num K_mod_num_approx] = kernel_m_singular_mod_numerical(zeta, m, CONSTS, p_m_vec, p_m_an_vec, n_beg_sum)

    a = CONSTS.a;
    k0 = CONSTS.k0;
    
    %p_m_vec = func_p_m(m, CONSTS, span, false);
    %n_beg_sum = n_begining_sum(CONSTS, p_m_vec);
    coef_m_mat = func_coef_m_of_eigenmodes(m, CONSTS, p_m_vec);
    
    sum = zeros(size(zeta), 1);
%     sum_an = zeros(size(zeta), 1);
       
    for i = n_beg_sum:size(p_m_vec, 1)

        p_n = p_m_vec(i);
        coef_n = coef_m_mat(:, i);
        norm = func_norm(m, CONSTS, coef_n, p_n);
        [~, ~, ~, ~, Ephi, ~] = func_HE_outside(1, m, CONSTS, coef_n, p_n); %func_HE_inside(a, m, CONSTS, coef_n, p_n);

        part_sum = zeros(size(zeta), 1);
        sum_n = zeros(size(zeta), 1);
        for ii = 1:size(zeta,1)
            zeta_i = zeta(ii);
            if(norm==0)
                part_sum(ii) = 0;
            else
                part_sum(ii) = (2*pi*a/norm)*(Ephi^2)*exp(-1i*k0*p_n*abs(zeta_i));
            end
            sum_n(ii) = sum_n(ii) + part_sum(ii);
        end
        
        sum = sum + sum_n;
        
    end
    
    K_mod_num = sum;
    
    clear i;
    
%     coef_m_an_mat = func_coef_m_of_eigenmodes(m, CONSTS, p_m_an_vec);
% 
%     for j = n_beg_sum:size(p_m_an_vec, 1)
%         
%         p_n_an = p_m_an_vec(j);
%         coef_n_an = coef_m_an_mat(:, j);
%         
%         norm_an = func_norm(m, CONSTS, coef_n_an, p_n_an);
%         [~, ~, ~, ~, Ephi_an, ~] = func_HE_outside(1, m, CONSTS, coef_n_an, p_n_an); %func_HE_inside(a, m, CONSTS, coef_n, p_n);
%         part_sum_an = zeros(size(zeta), 1);
%         sum_an_n = zeros(size(zeta), 1);
%         for jj = 1:size(zeta,1)
%             zeta_j = zeta(jj);
%             if(norm_an==0)
%                 part_sum_an(jj) = 0;
%             else
%                 part_sum_an(jj) = (2*pi*a/norm_an)*(Ephi_an^2)*exp(-1i*k0*p_n_an*abs(zeta_j));
%             end
%                sum_an_n(jj) = sum_an_n(jj) + part_sum_an(jj);
%         end
%         
%         sum_an = sum_an + sum_an_n;
%     end
    
    K_mod_num_approx = 0; %sum_an;
    
end

